563 

8.8A THE DISCRETE PROLATE SPHERIODAL FILTER AS A DIGITAL 
SIGNAL PROCESSING TOOL 

j. D. Mathews*, J. K. Breakall**, and G. K. Karawas* 

^Department of Electrical Engineering and Applied Physics, 

Case Institute of Technology, Case Western Reserve University, 

Cleveland, OH 44106 

**Lawreace Livermore National Laboratories, P.O.B. 808, Livermore, CA 94550 
ABSTRACT 

The discrete prolate spheriodal (DPS) filter is one of the class of non- 
recursive finite impulse response (FIR) filters. The DPS filter, first 
introduced by TUFTS and FRANCIS (1970), is superior to other filters in this 
class in that it has maximum energy concentration in the frequency passband and 
minimum ringing in the time domain. We give a mathematical development of the 
DPS filter properties, provide information required to construct the filter, and 
compare the properties of this filter with those of the more commonly used 
filters of the same class. We note that use of the DPS filter allows for 
particularly useful statements of data time/frequency resolution "cell” values 
and that overall it forms an especially useful tool for digital signal 
processing. 

INTRODUCTION 

The reduction of digitized data to final form usually involves what is 
broadly described as digital filtering. This process may be applied in either 
the domain in which sampling occurs (e.g. , time or space) or in the 
corresponding transform domain (e.g., frequency or spatial frequency). The net 
result of digital filtering is always to "smooth” or restrict the frequency 
content of the data. However, it is often very important that this smoothing or 
filtering process be done with minimum and well-understood distortion to the 
results in both domains. An example of an often misused filter is the equal 
weight ("square”) running average which "smooths" in the domain of application 
but results in a (sin x/x) ^ like passband shape in the transform domain. The 
(sin x/x) 2 behavior may have unacceptable sidelobe levels and/or spacing. 

The difficulty of minimizing distortion due to filtering starts with a 
statement of the criterion for judging when a filter is most satisfactory. We 
will follow the approach of TUFTS and FRANCIS (1970), PAPOULIS and BERTRAN 
(1972), and SLEPIAN (1978) in which an optimum finite (length) impulse response 
(FIR) non-recursive filter is derived on the basis of maximum energy concen- 
tration in the passband of the filter relative to energy in the total bandpass 
up to the Nyquist frequency o) N (=2tt/2t with t the sampling interval). Choice of 
this energy concentration criterion leads to filters described in terms of 
discrete prolate spheriodal sequences and wave functions. These sequences and 
wave functions have many desirable properties which were first investigated in 
the continuous case by SLEPIAN and POLLAK (1961), LANDAU and POLLAK (1961, 

1962), and SLEPIAN (1964). 

This paper is basically tutorial in that our major goals are to demonstrate 
the usefulness of FIR filters based on maximum passband energy concentration and 
to supply the information necessary for the "construction" of these filters. We 
do however introduce several apparently new results which prove useful in 
numerically solving the matrix equations which describe the filter. In the 
following section we develop the general equations for symmetric FIR digital 
filters and then obtain solutions based on the energy concentration criterion. 
Next we demonstrate that these filters have features which are desirable 
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especially when compared with more traditional filters, and finally, we present 
a summary and conclusions, 

MATHEMATICAL DERIVATIONS 

(a) A General Non-Recursive FIR Filter 

Using the time and frequency domains, a tapped delay line implementation of 
a non-recursive FIR filter acting on analog signal x(t) is shown schematically 
in Figure 1, The resultant output signal has the form 

K 

y(t) = l a x(t-m) (1) 

n--K 


where the coefficients a n are real and where the time reference is, for later 
convenience, at the center of the delay line. Assuming that the Fourier trans- 
form of x(t) -*-*X((i>) exists we Fourier transform (1) and find 


Y(u>) = H(oi) • X(w) 

where 


H(w) = f a e" jn “ T 
n=-K 

is the 'Voltage" transfer function and 
transfer function is 
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S(w) * H(oj) • H*(oj) = £ 
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IN — delay line- 



X( t+Kr) 


y(t) 


Figure 1. Block diagram of the non- recursive FIR filter 
described by Equation 1, 
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while the inverse Fourier transform of (3) is 


K 

h(t) * l a fi(t-m ) (5) 

n— K 


the filter impulse response with 6(t*) the unit impulse. Note that h(t) con- 
volved with x(t) yields equation (1). 

The filter described by (3) or (5) becomes a digital filter if we uniformly 
and instantaneously sample y(t) at intervals of t time. Then (1) becomes 


y j 


K 


l 

n=-K 


a x. 
n j-n 


( 6 ) 


where the index j refers to consecutive members of the set of sampled signals. 
If we restrict x(t) to be bandlimited to the Nyquist frequency (2tt/2t) or less, 
then the sampling/filtering process occurs without aliasing. 


(b) Choice of Coefficients 


Equations (3) or (5) describe the effects of the filter on an input signal 
given a particular set of coefficients {a }, These coefficients are often 
chosen such that the digital filter characteristics are similar to one of 
various common analog filters (e.g., the "ideal” filter). The process of 
synthesizing these filter characteristics often involves smoothing (windowing) 
of the resultant coefficient sequence to suppress ringing (BLACKMAN and TUKEY, 
1958; HAMMING, 1977), 


As mentioned in the introduction, we propose to choose the sequence {a } 
by maximizing the energy concentration of S(co) (equation (4)) in some interval 
[-w c ,co c J on [~w n ,od n ] where o) c is the filter "cutoff" frequency. Thus let 

w 0) 

a « / S(w) d iii/J S(w) dw (7) 

-a) ' -a) 

c n 


and find {a n } such that a is maximized. Substituting equation (4) for S(w) in 
(7) and performing the integrations yield 


a 


K 
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n ,m=~K 
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(n-m) tt 
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( 8 ) 


where e 


w /a) . If we define (a } as elements of vector a then (8) becomes 

c n n — 


T T 

a • E^. • a - aa • a = 0 

where each a. has (2K +1) elements and the elements of matrix JL are 


(9) 


e 

mn 


sin (n-m)Tre 
(n-m)7r 


(10) 


-K < m, n < K 

We are thus searching for a such that a is maximum. From Lagrange 
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multiplier theory, it is clear that (9) is satisfied by the set of (2K + 1) 
orthonormal eigenvectors .aW and the corresponding eigenvalues of the 
system 


(E 


W 


.a) 


(id 


,(°) 


where I is the identity matrix. We choose the elements of aA u/ the 
eigenvector associated with X^, the largest eigenvalue, as coefficients of our 
"optimum” filter. That is, aVo) maximizes equation (8) with a =* X^ for the 
K and e (i.e. , 


u) c > values in question. 


SLEPIAN (1978) refers to these eigenvectors and corresponding HCuO’s as 
discrete prolate spheriodal sequences and wave functions, respectively. He 
demonstrates that if we order the eigenvalues according to magnitude. 


1 ^ Xo ^ X]_ ^ ^ ^ ® 

and if a^*) is the n* th element of the eigenvector corresponding to \ ± then 


K 


l 

n=-K 


(i) a (j) 




1 i = j 
0 i 4 j 


Letting b be the n*th element of the eigenvector corresponding to eigen- 
value xy calculated for u> c ' * w N - u> c then (SLEPIAN, 1978; section 2.2) 


■ <-» i *» (i) 

(i) _ , t\ l n l , (2K-i) 


a_ 


(- 1 ) 1 


(12) 


|n| S K 


Slepian also notes that the discrete prolate spheroidal wave functions U(w) 
satisfy a Sturm-Liouville equation of the form 

[ co slot - cosire] ~ + t 2 [K(K+1)coso)t - 9] U = 0 (13) 

dco aw 

where 0 is the eigenvalue (different from X) corresponding to a particular 
wave function and all other parameters are as before. Substituting H(w) (=U) 
from (3) into (12) we find 

(1 K ' " e i ’ § (l) “ 0 (14) 

where for |n| < K 


e 1 


mn 


fl 

j (K+m) (K-m+1) 
- m^ cosve 

(K-m) (K+m+1) 
0 otherwise 


m » n - 1 > 0 
m - n 

m = n + 1 < K 


or equivalently 

~ (K+n) (K-n+l)a , + (n 2 cosire - 0 . )a + -~(K-n) (K+n+l)a = 0 
2 n-1 in 2 n+1 


(15) 


(16) 
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The matrix JE 1 is tridiagonal and thus equation (14) is much easier to 
solve numerically" than equation (11). Equation (14) yields the same normalized 
eigenvectors as (11) however the corresponding eigenvalues are different 
although we still seek the eigenvector corresponding to the largest eigenvalue 
9 q. Since we are interested only in a.(0) which is symmetric we can 
essentially halve the size of J2‘ by invoking symmetry (see equation (12)) 
which yields “ 

-0. a n + K(K+1 )a, = 0 
x 0 1 

plus equation (16) with 1 £ n < K. The compact form of (14) yields only eigen- 
values corresponding to "even 11 eigenvectors and the resultant .a*s may be 
normalized such that a T *a m 1. 

* 0, e 0 - K(K+1), and a n (|n| < K) * 
using the w c * 0 case and (12)) that 

a n = a n-l K + n “ 


From equation (14) we find that for 


while for 
, = and 


(i.e., e “ 1) we find £ 

(K - (n-1) ] 


1 < n < K 


We observe experimentally that K(K+1) > > K for all 0 < u < to„ and that 


( 0 ) ( 0 ) ( 0 ) 

a 0 > a l > a 2 > 


> a K (0) > 0 


1 > a. (0) (w ) > a. (0) (u 
-1 C - 1 c 


V 


where a. (oj c ) refers to the i'th component of the eigenvector corresponding 
to the largest eigenvalue for 0 < w c < cd™, The range constraints on the value 
of 0 q may be used in finding 0 q from (14). 

COMPARISON OF FILTER PROPERTIES 

In the previous section we showed that FIR filters with discrete prolate 
spheriodal sequences as coefficients (see equation (1)) have maximum energy 
concentration in the passband. These filters are known as discrete prolate 
spheriodal (DPS) filters and in this section we compare their properties with 
those of more traditional FIR filters. 

These "traditional" filters are all based on various windowing (weighting) 
functions applied to the coefficients of the infinite Fourier series representa- 
tion of the ideal low-pass filter. The coefficients of this series are 


C = sin mre 
n wn 


n = 0,±1,±2, ... 


e - to /to 
c n 


(17) 


and the windowing function W(n) is such that W(n) * 0 for |n| £ K thus 
truncating the series so that the resultant transfer function becomes 


H (to) - l W(n)C e _jn “ T 
W n=-K 


(18) 
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The W(n) (-K <_n <_K) are chosen to minimize in some sense the ringing (Gibb's 
phenomena) which results from truncation of the series representation of the 
ideal lowpass filter (HAMMING, 1977; Chap, 5* or RABINER and GOLD, 1975; Chap. 
3), In this discussion, we compare the DPS filter with the Hy(w) (equation 
(18)) of same length using rectangular, Kaiser, Von Hann, Dolph-Chebyshev, and 
Hamming windows. The W(n) corresponding to these windows are given or 
referenced in Table 1 and are discussed in more detail by HAMMING (1977) and 
RABINER et al. (1979). 

Figures 2 and 3 compare the frequency response (in dB) plotted versus 
normalized frequency f (Nyquist frequency f^ = 0,5) and impulse response 
envelope plotted versus time in samples for all six filters. For this compari- 
son we choose K « 15 (total length 31) and adjust filter parameters such that 
the -3 dB level of each filter occurs near f « 0.05. The filters appear in 
order of, in our opinion, acceptability with the rectangular windowed filter 
least acceptable and DPS filter most acceptable. 

The frequency domain distortion of a signal spectrum is just the filter 
frequency response as shown in Figures 2 and 3. However, the time domain 

distortions due to the impulse responses are not so obvious. In Figure 4, we 

show the effects of these filters in the time domain on a square wave plus 
Gaussian distributed noise. As expected, all filters except the DPS filter 
"ring 1 '. The DPS filter "smooths” the input signal while retaining a "smooth" 
frequency domain behavior. 

To aid in understanding the "evolution" of the DPS filter, we show in 
Figure 5 the DPS filter frequency response characteristics for normalized cut- 
off frequency f Q ® 0.2 as K is increased from 1 (3 coefficients) to 8 (17 

coefficients). Recall that in each case energy concentration in the passband 

(f < 0.2) is the maximum possible for the number of coefficients used. 

SUMMARY AND CONCLUSION 

We have investigated the properties of the discrete prolate spheroidal 
filter relative to other common FIR filters. To accomplish this comparison we 
first derive the general properties of FIR filters and then specialize these 

Table 1. Weighting coefficients of various windows 
(W(jn| > K) = 0 for all windows 



and we use a “ 3.395 (40 dB minimum stopband attenuation) 

1 + cos irn/K. 

VON HANN WINDOW : W(n) = 2 

DOLPH-CHEBYSHEV WINDOW: See HELMS (1968) and RABINER at al. (1979) 


HAMMING WINDOW ; W(n) = 0.54 + 0.46 cos (2 to/(K“1)') 
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properties to the case of maximum energy concentration in the filter passband. 
The condition of maximum energy concentration gives rise to an eigenvalue 
problem which has the discrete prolate spheroidal sequences (SLEPIAN, 1978) as 
eigenvector solutions* The eigenvector corresponding to the largest eigenvalue 
is used as the coefficient sequence in the DPS filter. 

SLEPIAN (1978) also points out that the prolate spheroidal wave functions 
satisfy a form of the Sturm-Liouville equation* This fact leads to a second 
eigenvalue/eigenvector problem which is much simpler to solve numerically than 
the first problem. We note that the elements of the eigenvector (a(°)) 
corresponding to the largest eigenvalue (0^ or X^) are symmetric around the 
center value (a Q ) and that 

a - (Q) 2a n <0) > 0 (1 < n < K). 
n-1 n 




Figure 2. Filter frequency response (dB) plotted versus normalized 
frequency (Nyquist frequency f N +) . 5 ) and corresponding impulse 
response envelope plotted versus time in samples for rectangular, 
Kaiser and von Hann windows of length 31 (see Equation 18 and Table 
1). In all cases the normalized 3 dB frequency is about 0.05. 
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When used, DPS filters have several features which we feel make them 
superior to more standard FIR filters. Specifically, the DPS filter is "smooth 
in both domains and thus distortion of the signal is minimal and easily 
described in both domains. Thus for example, after removing high frequency 
noise, a time domain feature may be more easily "restored 11 by deconvolving the 
DPS impulse response than by deconvolving the other sin x/x like impulse 
responses. Another useful feature of DPS filters is that descriptions of time 
and frequency resolutions of published data can be more readily appreciated and 
used by the reader* In fact, we suggest that the DPS filter might form a basis 
for standardization of information concerning the "resolution cells" of data. 
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Figure 3. Similar to Figure 2, but for Dolph-Chebysev and Hamming 
windows and the prolate spheroidal filter. The normalized cutoff 
frequency of the prolate spheroidal filter is 0,2. The filters 
in Figures 2 and 3 are ordered according to overall quality with, 
in our opinion, the rectangular window filter worst and prolate 
spheroidal filter best. 




effects of the various filters on a 
td noise, Note that all filters ex- 
'ring 1 ' and that the rectangular window 
window next while the Dolph-Chebysev, 
.milar in this case. 


Figure 5, The evolution of the prolai 
spheroidal filter frequency response 
as K increases from 1 through 8 for 
a normalized cutoff frequency of 0.: 
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